HACI-Fälle 2013, NC-Komposit des CI-Zeitpunktes¶
Für eine subjektive Einteilung der Fälle, ist es notwendig sich die CI-Zeipunkte anzusehen.
In [4]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
import xarray as xr
import glob
from standard_config import *
import sys
sys.path.append("{}/utils".format(local_home_path))
import load_satellite_data as lsd
import fixed_colourbar as fc
from analysis_tools import grid_and_interpolation as gi
import time
from IPython.display import clear_output
from l15_msevi import msevi as msv
import MSGtools as mst
from tqdm import tqdm
In [2]:
# nc_channels = ['VIS006','VIS008','IR_016']
channels = ['VIS006','VIS008','IR_016','IR_039','WV_062','WV_073','IR_087','IR_108','IR_120']
In [3]:
def cutout_box(data, row, col, box_size):
cutout = gi.cutout_field4box(data,(row,col),box_size)
return cutout
def normalise2range(data,dmin = None, dmax = None, new_min = 0, new_max = 1,gamma=1):
data = data**(1/gamma)
if dmin == None and dmax == None:
data_std = ((data - np.min(data)) / (np.max(data) - np.min(data)))
else:
data_std = ((data - dmin) / (dmax - dmin))
data_std = ((data - np.min(data)) / (np.max(data) - np.min(data)))
data_scaled = data_std * (new_max - new_min) + new_min
return data_scaled
def array_to_255(d, gamma = 1.0):
dat = d.copy()
# remove outliners for dat>1. or dat<0.
dat[dat>1.] = 1.
dat[dat<0.] = 0.
if not gamma == 1.0:
dat = dat ** (1./gamma)
# transform (0,1)-data to (0,255)-image-scale
rgb = 255 * dat
return rgb.astype('uint8')
def scale_array_min_max(array_data,range_min=0,range_max=1):
"""
Scales a array into the chosen range.
Inputs:
-------
array_data: numpy array of floats or integers, 2d
array to scale
range_min: int or float, default = 0
minimum value of the range to scale array to,
range_max: int or float, default = 1
maximum value of the range to scale array to,
Returns:
--------
scaled_array: numpy array of floats, 2d
"""
# get array extrema
array_min = np.min(array_data)
array_max = np.max(array_data)
# derive conversion parameters
a = (range_max - range_min) / (array_max - array_min)
b = range_max - a * array_max
# scale array
scaled_array = a * array_data + b
return scaled_array
def add_hrv_texture2nc(nc,hrv):
nc_lab = color.rgb2lab(nc)
l_min = nc_lab[...,0].min()
l_max = nc_lab[...,0].max()
l_hrv_scaled = nc_lab[...,0] * hrv
l_hrv_scaled = scale_array_min_max(l_hrv_scaled,l_min,l_max)
nc_lab[...,0] = l_hrv_scaled
return color.lab2rgb(nc_lab)
def gamma2(values,gamma):
vmin = np.min(values)
vmax = np.max(values)
vmean = (vmin + vmax)/2
below_mean = np.where(values < vmean)
above_mean = np.where(values >= vmean)
out = np.zeros(values.shape)
out[below_mean] = 128 - 128 * ((vmean - values[below_mean]) / (vmean - vmin))**(1/gamma)
out[above_mean] = 128 + 128 * ((values[above_mean] - vmean) / (vmax - vmean))**(1/gamma)
return out
def sat_rgb(red_channel,green_channel,blue_channel):
red = array_to_255(red_channel,1.0)
green = array_to_255(green_channel,1.0)
blue = array_to_255(blue_channel,1.0)
return np.dstack([red,green,blue]).astype("uint8")
def day_natural_composite(vis006_data,vis008_data,nir016_data):
(rmin, rmax, rgamma) = (0,1,1)
(gmin, gmax, ggamma) = (0,1,1)
(bmin, bmax, bgamma) = (0,1,1)
comp = sat_rgb(np.clip(nir016_data,rmin,rmax),
np.clip(vis008_data,gmin,gmax),
np.clip(vis006_data,bmin,bmax))
return comp
def day_microphysics_rgb(vis008_data,ir039_data,ir108_data):
(rmin, rmax, rgamma) = (0,1,1)
(gmin, gmax, ggamma) = (0, 0.6, 2.5)
(bmin, bmax, bgamma) = (203,323,1)
comp = sat_rgb(normalise2range(vis008_data, rmin, rmax, 0,1, rgamma),
normalise2range(ir039_data, gmin, gmax, 0,1, ggamma),
normalise2range(ir108_data, bmin, bmax, 0,1, bgamma))
return comp
def day_solar_rgb(vis008_data,ir016_data,ir039_data):
(rmin, rmax, rgamma) = (0,1,1.7)
(gmin, gmax, ggamma) = (0,0.7,1)
(bmin, bmax, bgamma) = (0,0.6,2.5)
comp = sat_rgb(normalise2range(vis008_data,rmin,rmax,0,1,rgamma),
normalise2range(ir016_data,gmin,gmax,0,1,ggamma),
normalise2range(ir039_data,bmin,bmax,0,1,bgamma))
return comp
def conv_rgb(vis006_data,ir016_data,ir039_data,wv062_data,wv073_data,ir108_data):
(rmin, rmax, rgamma) = (-30,0,1)
(gmin, gmax, ggamma) = (0,55,0.5)
(bmin, bmax, bgamma) = (-0.7,0.2,1)
r = normalise2range(wv062_data - wv073_data, rmin, rmax,0,1,rgamma)
g = normalise2range(ir039_data - ir108_data, gmin, gmax,0,1,ggamma)
b = normalise2range(ir016_data - vis006_data, bmin, bmax,0,1bgamma)
comp = sat_rgb(r,g,b)
return comp
def night_conv_rgb(ir039_data,ir108_data,ir120_data):
(rmin, rmax, rgamma) = (-4,2,1)
(gmin, gmax, ggamma) = (0,6,2)
(bmin, bmax, bgamma) = (243,293,1)
r = normalise2range(ir120_data - ir108_data, rmin, rmax, 0, 1, rgamma)
g = normalise2range(ir108_data - ir039_data, gmin, gmax, 0, 1, ggamma)
b = normalise2range(ir108_data, bmin, bmax, 0, 1, bgamma)
comp = sat_rgb(r,g,b)
return comp
def determine_IR039_reflectance(day,scan_type='rss'):
import datetime as dt
from l15_msevi import msevi as msv
from l15_msevi.msevi_config import _calibration_constants
import MSGtools as mst
s = msv.MSevi(day,chan_list=['IR_039','IR_108','IR_134'],scan_type=scan_type)
s.rad2bt()
# CO2-korrigierte Radianz
R_corr = (s.bt['IR_108'] - 0.25* (s.bt['IR_108'] - s.bt['IR_134'])**4) / (s.bt['IR_108']**4)
# CO2-korrigierte Helligkeitstemperatur
bt_039_korr = (s.bt['IR_039']**4 + R_corr)**0.25
# Konstanten
c1 = 1.19104e-5
c2 = 1.43877
nu_c = _calibration_constants[s.sat_type]['nu_c']['IR_039']
A = _calibration_constants[s.sat_type]['A']['IR_039']
B = _calibration_constants[s.sat_type]['B']['IR_039']
# Rückumrechnung von BT(10,8) zu L(10,8) mit Konstanten des IR039-Kanals
R_108 = (c1 * nu_c**3) / (np.exp((c2*nu_c) / (A*s.bt['IR_108'] + B)) - 1 )
# korrigierte thermische Komponente der Radianz
R_therm = R_108 * R_corr
# CO2-korrigierte Solarkonstante an der Atmosphärenoberseite
jday = int(day.strftime('%j'))
esd = 1. - 0.0167 * np.cos( 2*np.pi*(jday-3.) / 365 )
cos_sun_zenith = mst.get_cos_zen(day,scan_type=scan_type)
cos_sat_zenith = np.cos(np.deg2rad(mst.get_msg_sat_zen(day,sat_type=scan_type)))
solar_constant = 4.92 / esd**2
CO2_att_cloud_sat = np.exp(-(1 - s.rad['IR_039']))
CO2_att_cloud_sun = np.exp(-(1 - s.rad['IR_039']) * (cos_sun_zenith / cos_sat_zenith))
R_toa = solar_constant * cos_sun_zenith * CO2_att_cloud_sat * CO2_att_cloud_sun
# Reflektanz des IR039-Kanals
ref_039 = ((s.rad['IR_039'] - R_therm) / (R_toa - R_therm))
return bt_039_korr, ref_039
HACI-Falldaten laden und das Jahr 2013 auswählen¶
In [6]:
haci_objects = pd.read_csv("{}/HACI_bbox_data/haci-2008-2017-bbox-filtered.csv".format(local_data_path))
In [7]:
year = [int(pd.Timestamp(case.date).strftime("%Y")) for i, case in haci_objects.iterrows()]
In [9]:
haci_objects = haci_objects.assign(year=year)
In [10]:
haci_objects_2013 = haci_objects[haci_objects.year == 2013]
In [11]:
len(haci_objects_2013.index)
Out[11]:
Im Ganzen haben wir
NameError: name 'haci_objects_2013' is not defined
Fälle, die wir uns im Folgenden ansehen. Das sollte ausreichen, um einen repräsentativen Überblick über mögliche Fallklassen zu erhalten.Fälle durchlaufen und Ausschnitt des NC-Komposits für den CI-Zeitpunkt erstellen¶
In [12]:
pic_target = "{}/proj/2019-01_trackingstudie/pics/cutout_t0".format(local_home_path)
In [13]:
plt.switch_backend("Agg")
In [22]:
for i, case in tqdm(haci_objects_2013.iterrows(),total=len(haci_objects_2013.dt.values)):
#for i, case in haci_objects.iterrows():
try:
ci_time = pd.Timestamp("{}T{}".format(case.date,case.time)).to_pydatetime()
case_id = "{}_{}".format(ci_time.strftime("%Y%m%d"),case.id)
sat_data = lsd.load_satellite_data_multichannel(ci_time,channels)
cos_sun_zenith = mst.get_cos_zen(ci_time,scan_type='rss')
nc = day_natural_composite(cutout_box(sat_data['VIS006'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['VIS008'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_016'],case.msg_eu_l0,case.msg_eu_c0,51))
bt039, ref039 = determine_IR039_reflectance(ci_time)
mp = day_microphysics_rgb(cutout_box(sat_data['VIS008'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(ref039,case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_108'],case.msg_eu_l0,case.msg_eu_c0,51))
cv = conv_rgb(cutout_box(sat_data['VIS006'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_016'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_039'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['WV_062'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['WV_073'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_108'],case.msg_eu_l0,case.msg_eu_c0,51))
ncv = night_conv_rgb(cutout_box(sat_data['IR_039'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_108'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_120'],case.msg_eu_l0,case.msg_eu_c0,51))
#nc_cutouts[case_id] = nc
fig,ax = plt.subplots(1,4,figsize=(16,4))
ax[0].imshow(nc)
ax[0].axis("off")
ax[0].set_title("NC-RGB")
ax[1].imshow(mp)
ax[1].axis("off")
ax[1].set_title("Mikrophysik-RGB")
ax[2].imshow(cv)
ax[2].axis("off")
ax[2].set_title("Konvektions-RGB")
ax[3].imshow(ncv)
ax[3].axis("off")
ax[3].set_title("Nachtkonvektions-RGB")
plt.suptitle(case_id)
#plt.tight_layout()
plt.savefig("{}/{}.png".format(pic_target,case_id))
plt.close("all")
except Exception as e:
print("FEHLER {} bei Fall {}.".format(e,case_id))
continue
In [24]:
plt.switch_backend("Qt5Agg")
In [ ]:
In [23]:
sat_data
Out[23]:
In [25]:
ci_time = pd.Timestamp("{}T{}".format(case.date,case.time)).to_pydatetime()
case_id = "{}_{}".format(ci_time.strftime("%Y%m%d"),case.id)
In [26]:
ci_time
Out[26]:
In [40]:
i = 400
case = haci_objects_2013.iloc[i]
In [41]:
case
Out[41]:
In [54]:
import datetime as dt
In [62]:
# ci_time = pd.Timestamp("{}T{}".format(case.date,case.time)).to_pydatetime()
# case_id = "{}_{}".format(ci_time.strftime("%Y%m%d"),case.id)
ci_time = dt.datetime(2015,6,2,12,0)
sat_data = lsd.load_satellite_data_multichannel(ci_time,channels)
cos_sun_zenith = mst.get_cos_zen(ci_time,scan_type='rss')
In [56]:
sat_data['VIS006'] = sat_data['VIS006'] / cos_sun_zenith
sat_data['VIS008'] = sat_data['VIS008'] / cos_sun_zenith
sat_data['IR_016'] = sat_data['IR_016'] / cos_sun_zenith
In [59]:
%matplotlib inline
In [57]:
bt039, ref039 = determine_IR039_reflectance(ci_time)
nc = day_natural_composite(cutout_box(sat_data['VIS006'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['VIS008'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_016'],case.msg_eu_l0,case.msg_eu_c0,51))
mp = day_microphysics_rgb(cutout_box(sat_data['VIS008'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(ref039,case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_108'],case.msg_eu_l0,case.msg_eu_c0,51))
cv = conv_rgb(cutout_box(sat_data['VIS006'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_016'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_039'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['WV_062'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['WV_073'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_108'],case.msg_eu_l0,case.msg_eu_c0,51))
ncv = night_conv_rgb(cutout_box(sat_data['IR_039'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_108'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_120'],case.msg_eu_l0,case.msg_eu_c0,51))
ds = day_solar_rgb(cutout_box(sat_data['VIS008'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(sat_data['IR_016'],case.msg_eu_l0,case.msg_eu_c0,51),
cutout_box(ref039,case.msg_eu_l0,case.msg_eu_c0,51))
In [60]:
fig,ax = plt.subplots(1,5,figsize=(20,4))
ax[0].imshow(nc)
ax[0].axis("off")
ax[0].set_title("NC-RGB")
ax[1].imshow(mp)
ax[1].axis("off")
ax[1].set_title("Mikrophysik-RGB")
ax[2].imshow(cv)
ax[2].axis("off")
ax[2].set_title("Konvektions-RGB")
ax[3].imshow(ncv)
ax[3].axis("off")
ax[3].set_title("Nachtmikrophysik-RGB")
ax[4].imshow(ds)
ax[4].axis("off")
ax[4].set_title("Tag-Solar-RGB")
plt.suptitle(case_id)
Out[60]:
In [63]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
ax.imshow(conv_rgb(sat_data['VIS006'],
sat_data['IR_016'],
sat_data['IR_039'],
sat_data['WV_062'],
sat_data['WV_073'],
sat_data['IR_108']))
Out[63]:
In [ ]:
fig,ax = plt.subplots(1,1,figsize=(16,10))
ax.imshow(day_micro(sat_data['VIS006'],
sat_data['IR_016'],
sat_data['IR_039'],
sat_data['WV_062'],
sat_data['WV_073'],
sat_data['IR_108']))
In [64]:
bt039, ref039 = determine_IR039_reflectance(ci_time)
fig,ax = plt.subplots(1,1,figsize=(16,10))
ax.imshow(day_microphysics_rgb(sat_data['VIS008'],ref039,sat_data['IR_108']))
Out[64]:
In [631]:
from skimage import color
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
In [632]:
cv_rgb = conv_rgb(sat_data['VIS006'],
sat_data['IR_016'],
sat_data['IR_039'],
sat_data['WV_062'],
sat_data['WV_073'],
sat_data['IR_108'])
cv_lab = color.rgb2lab(cv_rgb)
In [633]:
mp_rgb = day_microphysics_rgb(sat_data['VIS008'],
ref_039,
sat_data['IR_108'])
mp_lab = color.rgb2lab(mp_rgb)
In [634]:
BuYl = LinearSegmentedColormap.from_list('BuYl', [[0,0,1,1],
[0.25,0.25,0.75,1],
[0.5,0.5,0.5,1],
[0.75,0.75,0.25,1],
[1,1,0,1]])
GnRd = LinearSegmentedColormap.from_list('GnRd', [[0,1,0,1],
[0.25,0.75,0,1],
[0.5,0.5,0,1],
[0.75,0.25,0,1],
[1,0,0,1]])
In [695]:
fig,ax = plt.subplots(2,2,figsize=(16,16))
axs= ax.ravel()
axs[0].imshow(cv_rgb)
axs[0].set_title("Convection RGB")
l_plot = axs[1].imshow(cv_lab[:,:,0],vmin=0,vmax=100,cmap='gray')
axs[1].set_title("L")
fc.colourbar(l_plot)
a_plot = axs[2].imshow(cv_lab[:,:,1],vmin=-70,vmax=70,cmap=GnRd)
fc.colourbar(a_plot)
axs[2].set_title("a")
b_plot = axs[3].imshow(cv_lab[:,:,2],vmin=-100,vmax=100,cmap=BuYl)
fc.colourbar(b_plot)
axs[3].set_title("b")
Out[695]:
In [646]:
fig,ax = plt.subplots(2,2,figsize=(16,16))
axs = ax.ravel()
axs[0].imshow(mp_rgb)
axs[0].set_title("Mikrophysik-RGB")
l_plot = axs[1].imshow(mp_lab[:,:,0],vmin=0,vmax=100,cmap='gray')
axs[1].set_title("L")
fc.colourbar(l_plot)
a_plot = axs[2].imshow(mp_lab[:,:,1],vmin=-100,vmax=100,cmap=GnRd)
fc.colourbar(a_plot)
axs[2].set_title("a")
b_plot = axs[3].imshow(mp_lab[:,:,2],vmin=-100,vmax=100,cmap=BuYl)
fc.colourbar(b_plot)
axs[3].set_title("b")
Out[646]:
